--------------------------------------------------------------------------------
-- This package is a simple implementation of common integer fractions.
--
-- Copyright (C) 2009 Vincent DIEMUNSCH.
-- GNU General Public License.
--
-- This program is free software: you can redistribute it and/or modify it under
-- the terms of the GNU General Public License as published by the Free Software
-- Foundation, either version 3 of the License, or (at your option) any later
-- version.
--
-- This file is distributed in the hope that it will be useful, but WITHOUT ANY
-- WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-- A PARTICULAR PURPOSE. See the GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License along with
-- this program.  If not, see <http://www.gnu.org/licenses/>.
--------------------------------------------------------------------------------

package body Rational_Numbers is

   function tdiv (left, Right : Integer) return Integer renames Standard."/";
   function tmod (left, Right : Integer) return Integer renames Standard."rem";
   -- Truncation Division, the standard Ada integer division.


   ----------------
   -- Simple GCD --
   ----------------

   function Simple_GCD (m, n : Natural) return Positive is
      a     : Natural := m;
      b     : Natural := n;
      New_a : Positive;
   begin
      while (b /= 0) loop
         New_a := b;
         b     := a mod b;       -- a,b = b, a mod b
         a     := New_a;
      end loop;
      return a;
   end Simple_GCD;

   -----------------
   -- Complex GCD --
   -----------------

   function Complex_GCD (m, n : Natural) return Positive is
      a        : Natural;
      b, New_a : Positive;
   begin

      if n = 0 then return m; end if;
      if m = 0 then return n; end if;

      a := m;
      b := n;

      while (a mod 2) = 0 and then (b mod 2) /= 0 loop a := tdiv(a,2); end loop;
      while (b mod 2) = 0 and then (a mod 2) /= 0 loop b := tdiv(b,2); end loop;

      if a /= 1 and then b /= 1 then

         -- initial reduction using a mod b, worthwhile if b is much smaller
         -- than a. Otherwise, it is faster to iterate 2 or 3 times (a-b)
         if a > 4*b then
            a := a mod b;                 -- a,b = a mod b, b
            if a = 0 then
               return b;                  -- GCD found !
            else
               while (a mod 2)  = 0 and then
                     (b mod 2) /= 0 loop  -- strip factors of 2 from a
                  a := tdiv(a,2);             -- (to keep it odd !)
               end loop;
               null;
            end if;
         end if;

         -- Main loop, see gmp manual for details, or Wikipedia !
         while a /= b loop
            -- a,b = abs(a-b), min(a,b) :
            New_a := abs(a-b);
            if a < b then b := a; end if;
            a := New_a;
            -- strip factors of 2 from a :
            while (a mod 2) = 0 and then (b mod 2) /= 0 loop
               a := tdiv(a,2);
            end loop;
         end loop;

         return a;
      else
         return 1;
      end if;
   end Complex_GCD;

   ------------------
   -- Canonicalize --
   ------------------

   procedure Canonicalize (F : in out Rational) is
      a : Positive;
   begin
      -- trivial cases :
      if    F.num  = 0 then F.den := 1;
      elsif F.num /= 1 then

         -- After the following  2 is never a common divisor anymore.
         -- we strongly hope that the compiler knows that in binary (a mod 2) is
         -- simply (a and 1) and that (a / 2) is "shift right a". Otherwise
         -- performances will be poor !
         while (F.num mod 2) = 0 and then (F.den mod 2) = 0 loop
            F.num := tdiv(F.num,2);
            F.den := tdiv(F.den,2);
         end loop;

         -- Find Greatest Common Divisor : a --
         a := Complex_GCD (abs F.num, F.den);
         Pragma Assert (a = Simple_GCD (abs F.num, F.den)); -- DEBUG ONLY !

         -- Simplify :
         if a > 1 then
            F.num := tdiv(F.num,a);
            F.den := tdiv(F.den,a);
         end if;
      end if;
   end Canonicalize;

   ---------
   -- "/" --
   ---------

   function "/" (Left, Right : integer) return Rational is
      Result : Rational;
   begin
      if Right >= 0 then Result := ( Left,  Right); -- Constraint_Error if Right = 0
                    else Result := (-Left, -Right);
      end if;
      Canonicalize (Result);
      return Result;
   end "/";

   ---------
   -- Num --
   ---------

   function Num (F : Rational) return Numerator is
   begin
      return F.Num;
   end Num;

   ---------
   -- Den --
   ---------

   function Den (F : Rational) return Denominator is
    begin
      return F.Den;
   end Den;

   ----------------------------
   -- Set_Canonical_Rational --
   ----------------------------

   function Set_Canonical_Rational (Left : Numerator; Right : Denominator) return Rational is
    begin
      return (Left, Right);
   end Set_Canonical_Rational;

   ---------
   -- GCD --
   ---------

   function GCD (m, n : Natural) return Positive renames Complex_GCD;

   -----------
   -- Image --
   -----------
   function Image (F : Rational) return String is
   begin
      return Integer'Image(F.num) & " /" & Integer'Image (F.den);
   end Image;


   ---------
   -- "=" --
   ---------

   function "=" (Left, Right : Rational) return Boolean is
   begin
      if Left.den = Right.den and then Left.num = Right.num then
         return True;
      else
         return (Left.num*Right.den = Right.num*Left.den);
      end if;
   end "=";

   ---------
   -- "<" --
   ---------

   function "<" (Left, Right : Rational) return Boolean is
   begin
      if Left.den = Right.den and then Left.num < Right.num then
         return True;
      else
         return (Left.num*Right.den < Right.num*Left.den);
      end if;
   end "<";

   ----------
   -- "<=" --
   ----------

   function "<=" (Left, Right : Rational) return Boolean is
   begin
      if Left.den = Right.den and then Left.num <= Right.num then
         return True;
      else
         return (Left.num*Right.den <= Right.num*Left.den);
      end if;
   end "<=";

   ---------
   -- ">" --
   ---------

   function ">" (Left, Right : Rational) return Boolean is
   begin
      if Left.den = Right.den and then Left.num > Right.num then
         return True;
      else
         return (Left.num*Right.den > Right.num*Left.den);
      end if;
   end ">";

   ----------
   -- ">=" --
   ----------

   function ">=" (Left, Right : Rational) return Boolean is
   begin
      if Left.den = Right.den and then Left.num >= Right.num then
         return True;
      else
         return (Left.num*Right.den >= Right.num*Left.den);
      end if;
   end ">=";

   ---------
   -- "+" --
   ---------

   function "+" (Right : Rational) return Rational is
   begin
      return Right;
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (Right : Rational) return Rational is
   begin
      return Rational'(-Right.num, Right.den);
   end "-";

   -----------
   -- "abs" --
   -----------

   function "abs" (Right : Rational) return Rational is
    begin
      return Rational'(abs(Right.num), Right.den);
   end "abs";

   ---------
   -- "+" --
   ---------

   function "+" (Left, Right : Rational) return Rational is
      Result : Rational;
   begin
      if Left.den = Right.den then
         Result.num := Left.num + Right.num;
         Result.den := Left.den;
      else
         Result.num := Left.num*Right.den + Right.num*Left.den;
         Result.den := Left.den * Right.den;
      end if;
      Canonicalize (Result);
      return Result;
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (Left, Right : Rational) return Rational is
      Result : Rational;
   begin
      if Left.den = Right.den then
         Result.num := Left.num - Right.num;
         Result.den := Left.den;
      else
         Result.num := Left.num*Right.den - Right.num*Left.den;
         Result.den := Left.den * Right.den;
      end if;
      Canonicalize (Result);
      return Result;
   end "-";

   ---------
   -- "*" --
   ---------

   function "*" (Left, Right : Rational) return Rational is
      Result : Rational := (Left.num * Right.num, Left.den * Right.den);
   begin
      Canonicalize (Result);
      return Result;
   end "*";

   ---------
   -- "/" --
   ---------

   function "/" (Left, Right : Rational) return Rational is
      Result : Rational;
   begin
      if Right.num = 0 then raise Constraint_Error; end if;
      Result := (Left.num * Right.den, Left.den * Right.num);
      Canonicalize (Result);
      return Result;
   end "/";

   ---------
   -- inv --
   ---------

   function inv (Right : Rational) return Rational is
   begin
      if Right.num = 0 then raise Constraint_Error; end if;
      return (Right.den, Right.num);
   end inv;

   ----------
   -- "**" --
   ----------

   function "**" (Left : Rational; Right : Integer'Base) return Rational is
   begin
      -- No need to simplify the result if the given Rational was simplified.
      return (Left.num ** Right, Left.den ** Right);
   end "**";

   ---------
   -- "+" --
   ---------

   function "+" (Left : Integer'Base; Right : Rational) return Rational is
      F : Rational := (Left * Right.den + Right.num, Right.den);
   begin
      Canonicalize (F);
      return F;
   end "+";

   function "+" (Left : Rational; Right : Integer'Base) return Rational is
   F : Rational := (Left.num + Right * Left.den, Left.den);
   begin
      Canonicalize (F);
      return F;
   end "+";

   ---------
   -- "-" --
   ---------

   function "-" (Left : Integer'Base; Right : Rational) return Rational is
      F : Rational := (Left * Right.den - Right.num, Right.den);
   begin
      Canonicalize (F);
      return F;
   end "-";

   function "-" (Left : Rational; Right : Integer'Base) return Rational is
   F : Rational := (Left.num - Right * Left.den, Left.den);
   begin
      Canonicalize (F);
      return F;
   end "-";


   ---------
   -- "*" --
   ---------

   function "*" (Left : Integer'Base; Right : Rational) return Rational is
      F : Rational := (Left, Right.den);
   begin
      Canonicalize (F); -- only the part that is new is simplified !
      F.num := F.num * Right.num;
      return F;
   end "*";

   ---------
   -- "*" --
   ---------

   function "*" (Left : Rational; Right : Integer'Base) return Rational is
      F : Rational := (Right, Left.den);
   begin
      Canonicalize (F); -- only the part that is new is simplified !
      F.num := F.num * Left.num;
      return F;
   end "*";

   ---------
   -- "/" --
   ---------

   function "/" (Left : Integer'Base; Right : Rational) return Rational is
      F : Rational;
   begin
      if    Right.num = 0 then raise Constraint_Error;
      elsif Right.num < 0 then F := (-Left, -Right.num);
                          else F := ( Left,  Right.num);
      end if;
      Canonicalize (F); -- only the part that is new is simplified !
      F.num := F.num * Right.den;
      return F;
   end "/";

   ---------
   -- "/" --
   ---------

   function "/" (Left : Rational; Right : Integer'Base) return Rational is
      F : Rational;
   begin
      if    Right = 0 then raise Constraint_Error;
      elsif Right < 0 then F := (-Left.num, -Right);
                      else F := ( Left.num,  Right);
      end if;
      Canonicalize (F); -- only the part that is new is simplified !
      F.den := F.den * Left.den;
      return F;
   end "/";

   ---------
   -- div --
   ---------

   function div (m, n : integer) return Integer is
      Result : Integer := tdiv(m,n);
   begin
      if tmod(m,n) < 0 then
         if n > 0 then Result := Result -1;
                  else Result := Result +1;
         end if;
      end if;
      return Result;
   end div;

   ---------
   -- mod --
   ---------

   function "mod" (m,n : integer) return Integer is
      Result : Integer := tmod(m,n);
   begin
      if Result < 0 then
         if n > 0  then Result := Result + n;
                   else Result := Result - n;
         end if;
      end if;
      return Result;
   end "mod";


   -----------
   -- Floor --
   -----------

   function Floor (F : Rational) return Integer is
     Result : Integer;
   begin
      if F.den = 1 then return F.num;
      else
         Result := tdiv(F.num,F.den);
         if (F.num < 0) /= (F.Den < 0) and then tmod(F.num,F.den) /= 0 then Result := Result - 1; end if;
         Pragma Assert (Result = Integer( Float'Floor (Float(F.num) / Float(F.den))) ); -- DEBUG ONLY !
         return Result;
      end if;
   end Floor;


   -------------
   -- Ceiling --
   -------------

   function Ceiling (F : Rational) return Integer is
     Result : Integer;
   begin
     if F.den = 1 then return F.num;
      else
         Result := tdiv(F.num,F.den);
         if Result >= 0 and then tmod(F.num, F.den) /= 0 then Result := Result + 1; end if;
	 Pragma Assert (Result = Integer( Float'Ceiling (Float(F.num) / Float(F.den))) ); -- DEBUG ONLY !
	 return Result;
	 end if;
   end Ceiling;


   --------------
   -- Rounding --
   --------------

   function Rounding (F : Rational) return Integer is
     Result : Integer;
   begin
      if F.den = 1
      then return F.num;
      else
         Result := tdiv(F.num,F.den);
         if tmod(F.num, F.den) /= 0 then
            Result := Result + 2*tdiv (tmod(F.num,F.den), F.den);
         end if;
         Pragma Assert (Result = Integer( Float'Rounding (Float(F.num) / Float(F.den))) ); -- DEBUG ONLY !
         return Result;
      end if;
   end Rounding;


   -----------------------
   -- Unbiased_Rounding --
   -----------------------

   function Unbiased_Rounding (F : Rational) return Integer is
     Result, Twice_Rem : Integer;
   begin
      if F.den = 1 then return F.num;
      else
         Result    := tdiv(F.num,F.den);
         Twice_Rem := 2 * tmod(F.num,F.den);
         if Twice_Rem /= 0 and then not (abs(Twice_Rem) = F.den and then tmod(Result,2) = 0) then
            Result := Result + tdiv(Twice_Rem,F.den);
         end if;
         Pragma Assert (Result = Integer( Float'Unbiased_Rounding (Float(F.num) / Float(F.den))) ); -- DEBUG ONLY !
         return Result;
      end if;
   end Unbiased_Rounding;


   ----------------
   -- Truncation --
   ----------------

   function Truncation (F : Rational) return Integer is
   begin
     return Standard."/"(F.num, F.den);
   end Truncation;


   --------------
   -- To_Float --
   --------------

   function To_Float (Q : Rational) return F is
   begin
      return F(Q.num) / F(Q.den);
   end To_Float;

   --------------------
   -- To_Fixed_Point --
   --------------------

   function To_Fixed_Point (Q : Rational) return F is
   begin
      return F(Q.Num) / F(Q.Den);
   end To_Fixed_Point;

end Rational_Numbers;
